##################### EXAMPLE - Raw simulation ############################
set.seed(9)         # to specify the initial seed
Nsim=10^4         # number of random variables
U=runif(Nsim)   #random numebrs from a uniform
X=-log(U) 
Y=rexp(Nsim)
par(mfrow=c(1,2))
hist(X,freq=F,main="Exp from Uniform")
hist(Y,freq=F,main="Exp from the R's function")
# We can also simulate a multivariate normal variable using univariate normals:
# 1) Cholesky decomposition of S = AA'
# 2) If Y is distributed as a N(0, I) then AY is distributed as a N(0, S)
# There is an R package that replicates those steps, called rmnorm.

##################### EXAMPLE - Solve the definite integral ############################
#the definite integral is x^3/3+3x^2/2+x+c
#therefore the definite integral between 0 and 1 is 1/3+3/2+1-0=17/6=2.833
#between 0 and 4 is 64/3+24+4-0=148/3=49.333

#MC estimation of integral between 0 and 1
TI=2.833
n=c(20,200,2000)
stima=0
for(i in 1:length(n)){
	x=runif(n[i])
	inte=0
	for(j in 1:length(x)){
		inte=inte+x[j]^2+3*x[j]+1
		}
	stima[i]=inte/n[i]
	}
stima
TI
#the higher the value of n, the closer the estimation is to the true integral

#MC estimation of integral between 0 and 4
TI=49.333
n=c(20,200,2000,20000,2000000)
stima=0
for(i in 1:length(n)){
	x=runif(n[i],0,4)
	inte=0
	for(j in 1:length(x)){
		inte=inte+x[j]^2+3*x[j]+1
		}
	stima[i]=4*(inte/n[i])
	}
stima
TI
#the higher the value of n, the closer the estimation is to the true integral